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ABSTRACT 

Continued radial velocity monitoring of the nearby M4V red dwarf star G J 876 
with Keck/HIRES has revealed the presence of a Uranus-mass fourth planetary 
companion in the system. The new planet has a mean period of Pg = 126.6 days 
(over the 12.6-year baseline of the radial velocity observations), and a minimum 
mass of nie sin = 12.9 ± 1.7 Af®. The detection of the new planet has been en- 
abled by significant improvements to our radial velocity data set for GJ 876. The 
data have been augmented by 36 new high-precision measurements taken over 
the past five years. In addition, the precision of all of the Doppler measurements 
have been significantly improved by the incorporation of a high signal-to-noise 
template spectrum for GJ 876 into the analysis pipeline. Implementation of the 
new template spectrum improves the internal RMS errors for the velocity mea- 
surements taken during 1998-2005 from 4.1 ms~^ to 2.5 ms~^. Self-consistent, 
N-body fits to the radial velocity data set show that the four-planet system has an 
invariable plane with an inclination relative to the plane of the sky of i = 59.5°. 
The fit is not significantly improved by the introduction of a mutual inclination 
between the planets "b" and "c," but the new data do confirm a non-zero eccen- 
tricity, = 0.207 ±0.055 for the innermost planet, "d." In our best-fit coplanar 
model, the mass of the new component is rUe — 14.6 ± 1.7 M^. Our best-fitting 
model places the new planet in a 3-body resonance with the previously known 
giant planets (which have mean periods of Pc = 30.4 and Pb = 61.1 days). The 
critical argument, (/'Laplace — ~ 3 A;, -|- 2Xg, for the Laplace resonance librates 
with an amplitude of Ac/^Lapiace = 40±13° about </?Lapiace = 0°. Numerical integra- 
tion indicates that the four-planet system is stable for at least a billion years (at 
least for the coplanar cases). This resonant configuration of three giant planets 
orbiting an M-dwarf primary differs from the well-known Laplace configuration 
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of the three inner Gahlean sateUites of Jupiter, which are executing very small 
librations about </5Lapiacc = 180°, and which never experience triple conjunctions. 
The GJ 876 system, by contrast, comes close to a triple conjunction between the 
outer three planets once per every orbit of the outer planet, "e." 

Subject headings: stars: GJ 876 - planetary systems - planets and satellites: general 
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Introduction 



The planetary system orbiting the nearby M4V star GJ 876 (HIP 113020) has proven 
to be perhaps the most remarkable entry in the emerging galactic planetary census. This 
otherwise unassuming red dwarf has produced the first example of a giant planet orbiting 
a low-mass star, the first instance of a mean-motion resonance among planets, the first 
clear-cut astrometric detection of an extrasolar planet, and one of the first examples of a 
planet in the hitherto unknown mass regime between Earth and Uranus. 

The radial velocity (RV) variations of GJ 876 have been monitored by us, the 
Lick-Carnegie Exoplanet Survey team (LCES), at the Ke ck I telescope using the High 



Resolution Echelle Spectrograph (HIRES) for 12.6 years. 



Bean &: Seifahrt 



Rivera et al. 



(120091) and 



(120051 ) (R05 henceforth) give detailed discussions of the discoveries and some 



of the studies made for this system. Here, we summarize some of that history. 



Marcy et al. 



(1l998h and 



Delfosse et al. 



(jl998[ ) announced the first companion, "b." 



They found that it has an orbital period (Pb) of ~ 61 days and a minimum mass {mi) sin it,) 
of ~ 2.1Mjup and that it produced a reflex barycentric velocity variation of its host star 
of amplitude Kh ~ 240 ms^^. After 2.5 more years of additional Doppler monitoring. 



Marcy et al. 



(120011 ) announced the discovery of a second giant planet in the system. 
This second companion, "c," has an orbital period, Pc ~ 30 days, nicsmic ~ 0.56 Mjup, 
and Kc ~ 81 ms~^, and upon discovery of t he sec ond planet, the parameters of the first 



planet were tangibly revised. 



Marcy et al. 



(|2001[ ) modeled the two-planet system with 
non-interacting Keplerian orbits, and noted both that there was room for improvement in 
the quality of the fit and that the system's stability depended on the initial positions of the 
planets. 

To a degree that has not yet been observed for any other planetary system, the mutual 
perturbations between the planets are dynamically significant over observable time scales. 
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This fortuitous situation arises because the minimum masses are relatively large compared 
to the star's mass, because the planetary periods are near the 2:1 commensurability, and 
because their orbits are confined to a small region around the star. 



Laughlin &: Chambers I (1200 ll ) and iRivera fc Lissauer I (1200 ll ) independently developed 
self-consistent "Newtonian" fitting schemes which incorporate the mutual perturbations 
among the planets in fitting the RV data. The inclusion of planet-planet interactions 
resulted in a substantially improved fit to the RV data, and suggested that the coplanar 



inclination o 



Nauenberg 



the b-c pair lies near i = 50°. These papers were followed by an article by 



(I2OO2I ) who described a similar dynamical fitting method, but, in contrast. 



found a system inclination near i = 90°. 



Soon after the announcement of planet "c," iBenedict et al. I (120021 ) used the Fine 
Guidance Sensor on the Hubble Space Telescope to detect the astrometric wobble induced by 
planet "b," which constituted the first unambiguous astrometric detection of an extrasolar 
planet. Their analysis suggested that the orbital inclina tion of plane t "b" i s close to edge-on 



84° ±6°), a result that agreed with the model 



in confiict with the results of 



Laughlin fc Chambers 



by 



Nauenberg 



(120011 ) and 



(120021). but whi c h was 



Rivera fc Lissauer 



(j200l|). 



R05 analyzed an updated RV data set that included new Doppler velocities obtained 
at Keck between 2001-2005. By adopting self-consistent fits, they were able to announce 
the discovery of a third companion, "d," with period Pd ~ 1.94 days, m^sinzd ~ 5.9 M^, 
and Kd ~ 6.5 ms~^. Using only the RV data, they were able to constrain the inclination 
of the system (assuming all three planets are coplanar) relative to the plane of the sky to 
50 ± 3°. For this inclination, the true mass of the third companion is ~ 7.5 M^. 



Bean fc Seifahrt 



(I2OO9I ) perforr ned self- consiste n t fitti ng of both the Keck RV data 



from R05 and the astrometry from 



Benedict et al. 



(I2OO2I ). They found the that the 



joint data set supports a system inclination i = 48.9]*;} go, which is in agreement with the 
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i = 50 ± 3° value published by R05. Bean & Seifahrt's analysis (see, e.g. their Figure 1) 
indicates that the best-fit xl is not significantly aff ected by the inc l usion of the astrometric 



Benedict et al. 



(I2OO2I ) is not in conflict 



data, and that the astrometric detection of "b" by 
with the significantly inclined system configuration. Bean & Seifahrt's work also resulted 
in the first determination of a mutual inclination between the orbits of two planets in an 
extrasolar planetary system. By allowing the inclinations, ib, and ic, and the nodes, Qb, and 
Qc to float as free parameters in their three-planet dynamical flt, they derived a mutual 



inclination, cos ^bc = cos ib cos ic + sin ib sin ic cos {Qc — ^b), of ^bc = 5.0 



+3.9° 

-2.3°- 



Correia et al. 



(I2OO9I ) also performed self-consistent, mutually inclined flts using the 
RVs from R05 plus 52 additional high precision RV measurements taken with the HARPS 
spectrograph. They conflrmed the presence of companio n "d," and further sh owed that it 



Correia et al. 



f l2009h determined 



has a signiflcant orbital eccentricity, = 0.139 ± 0.032. 
inclinations ib = 48.93° ± 0.97° and ib = 48.07° ± 2.06°. Their orbital fit yields a mutual 
inclination ^bc = 1.00°, consistent to wit hin measurement unc ertainty with a coplanar 
system configuration. We became aware of L 



Correia et al. 



( I2OO9I ) while preparing this article 



for publication. A. Correia (2009, personal communication) graciously shared the HARPS 
radial velocity data on which their paper is partly based. We have carried out a preliminary 
analysis that indicates that our new results are not in confiict with these HARPS data. The 
HARPS RVs do not add much to constrain the system's dynamics because of poor phase 
coverage at the period of the forth planet to be discussed below. Throughout this article, 
however, our analysis and results are based on the Keck RVs alone. 



In addition to R05, several other studies, such as 



West 



et al. (2005), and 



Shankland et al. 



(I1999I) 



West 



(1200l[ l. Laughlin 



( I2OO6I ). have examined the possibility of transits 



occurring in this system. The star has also been included in survey s in which it was prob ed 



for nearby companions or a circumstellar dust disk. These include 



Leinert et al. 



fll997h . 
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Patience et al. (2002) 



Hinz et al. I (120021) . iTrilling et al. I (120001) . iLestrade et al. I (120061) 



Shankland et al. 



(120081), and 



Bryden et al. 



(l2009h 



Luhman fc Jayawardhana 



(l2002h 



were able to use Keck A O to place limits o n 



separations of 1-10 AU. 



Seager fc Deming 



pote ntial brown dwarfs in the system at 



(120091 ) used the Spitzer spacecraft to search 



for thermal em ission from planet "d." The star was also included in the Gemini Deep 



Planet Survey (iLafreniere et al. 



20071 ). In addition to some of the work above, a variety of 



stu dies have examined t he dy n amics, and/or long-term stab i 



are 



Kinoshita &: Nakai 



(|20m 



and 



Snellgrove et al. 



2001) 



Haghighipour et al. 



(1200 ih 



Gozdziewski fc Macieiewski 



( 2OO3I ). 



Ji et al. 



(2002) 



ity of the system. Examples 



Jones et a. 



Beauge fc Michtchenko 



(2001) and 



mm 



Gozdziewski et al. 



Zhou &: Sun 



(|2003[ ) 



(120031 ) present an analytic 



model which does a v ery good job a t 



Murrav et al. 



(l2003h. 



(2002), 



Lee fc Peale 



Beauge et al. 



Lee &: Thommes 



fe003h 



repr o ducing the dynamical evolution of the system . 



(12002 ) 



Kley et al. 



. Chiang et a 




( 


2005 


), 


Thommes & Lissauer 


(2004 




Lee 


( 


2004 


), 


Kiev et al. 


( 


2005 


), and 



(I2OO9I ) examined scenarios of how the system may have come to be in 
its current dynamical configura tion. Even interior structure models of GJ 876 d have been 



introduced (IValencia et al. 



20071 ) 



In this work, we present a detailed examination of our significantly augmented and 
improved HIRES RV data set from the Keck observatory. Our analysis indicates the 
presence of a fourth, Uranus-mass planet in the GJ 876 system. The new planet is in 
a Laplace resonance with the giant planets "b" and "c," and the system marks the first 
example of a three-body resonance among extrasolar planets. We expect that the newly 
resolved presence of the fourth planet will provide additional constraints on the range of 
formation scenarios that could have given rise to the system, and will thus have a broad 
impact on current theories for planetary formation. The plan for this paper is as follows: In 
§ El we describe the new velocities and the methods used to obtain them. In § 131 we present 
our latest two- & three-planet coplanar fits. In § HI we thoroughly analyze the residuals of 
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our best three-planet coplanar fits. In § [5|, we discuss our four-planet fits. In § [6] we explore 
the possibility of constraining the libration amplitude of the critical angle for the Laplace 
resonance. In § [7] we explore the stability of additional planets in the system. Finally, in 
§ [8], we finish with a summary and a discussion of our findings. 



2. Radial Velocity Observations 



he st ellar characteris t ics of GJ 876 have been described previously i n iMarcv et al. 



fll998h and 



Laughlin et al. 



(120051 ). It has a Hipparcos distance of 4.69 pc (jPerryman et al. 



19971 ). Its luminosity is 0.0124 Lq. As in previous studies, we adopt a stellar mass of 0.32 M c^ 



Henry fc McCarthy 







and a radius of 0.3 i?© based on the mass-luminosity relationship of 
(119931 ). We do not incorporate uncertainties in the star's mass (0.32 ± 0.03 M©) into the 
uncertainties i n planetary masses and semi-major axes quoted herein. The age of the star 



exceeds 1 Gyr (IMarcy et al. 



19981) 



We searched for Doppler variability using repeated, high resolution sp ectra with 



resolving power R ^ 70000, obtained with Keck/HIRES (IVogt. et al. Ill994h . The Keck 



spectra span the wavelength range from 3700-8000 A. An iodine absorpt ion cell provides 



wave l ength calibration a nd the instrumental profile from 5000 to 6200 A ( Marcy Sz Butler 



1992 



Butler et al. 



19961 ). Typical signal-to-noise ratios are 100 per pixel for GJ 876. At 



Keck we routinely obtain Doppler precision of better than 3-5 ms~^ for V= 10 M dwarfs. 
Exposure times for GJ 876 and other V = 10 M dwarfs were typically 8 minutes. Over the 
past year, these typical exposure times have been raised to 10 minutes. 

The internal uncertainties in the velocities are judged from the velocity agreement 
among the approximately 700 2-A chunks of the echelle spectrum, each chunk yielding an 
independent Doppler shift. The internal velocity uncertainty of a given measurement is the 
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uncertainty in the mean of the ~ 700 velocities from one echelle spectrum. 

We present results of N-body fits to the RV data taken at the W. M. Keck telescope 
from 1997 June to 2010 January. The 162 measured RVs are listed in Table [H The 
median of th e internal uncer t aintie s is 2.0 ms~^. Comparison of these velocities with those 



presented in 



Laughlin et al. 



(120051 ) and in R05 shows significant changes (typically 3-10 
ms~^) in the velocities at several observing epochs. R05 discuss the primary reasons for the 
improvements — a CCD upgrade and an improved data reduction pipeline. Additionally, 
many of the velocities presented here are the result of averaging multiple exposures over 
two-hour time bins. More importantly, the present RV data set is based on a new "super" 
template. For the new template, ten exposures, most of which were < 900 sec, were 
obtained under very good observing conditions. In contrast, the template used to determine 
the RVs and uncertainties in R05 was based on only three 900-sec exposures taken under 
relatively poor observing conditions. As a result, the new template has brought the median 
internal uncertainty down to the current 2.0 ms~^ from the value of 4.1 ms~^ in R05. 
Additionally, the new template has resulted in significant changes in the old published 
RVs. Note that since the last observation listed in R05, 36 more high quality observations 
have been obtained for this system since 2004 December. Our LCES survey assigns high 
observing priority to GJ 876, given the potential for further characterizing the intricate 
dynamical configuration between the 30- and 61-day planets and because of the possibility 
of detecting additional planets in the system. 



3. Two- & Three-Planet Coplanar Fits 

We used similar procedures to those described in R05 to produce updated two- and 
three-planet dynamical fits to the GJ 876 Doppler velocities. The epoch for all the fits in 
this work is JD 2450602.093. As in R05, we initially held the eccentricity of companion 
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"d" (crf) fixed. We find, however, that fitting for results in a significant improvement 
in all the fits which include companion "d." We therefore allow to fioat in our final 
fits. Uncertainties in the tables listing our best-fit jacobi parameters are based on 1000 
bootstrap realizations of the RV data set. We take the uncertainties to be the standard 
deviations of the fitted parameters to the 1000 tria ls. The fitting algo rithm is a union of 



Levenberg-Marquar dt (LM) minimiza tion (§15.5 of IPress et al. I (119921 )) and the Mercury 



integration package (jChambers 



19991 ). 



Our dynamical fits can be integrated to determine the libration amplitudes of the 
critical arguments of the 2:1 MMR between "b" and "c," and to monitor the linear secular 
coupling between the pericenter longitudes. For each fit, we integrated the system for 
300,000 days (~821 years) with a timestep of 0.05 day. For all simulation s in this work, we 



used the Hybri d symplectic algorithm in the Mercury integration package (IChambers 



modified as in 



Lissauer fc Rivera 



19991 ) 



( 120011 ) to include the partial first order post-Newtonian 
correction to the central star's gravitational potential. The three relevant angles are 
V^cfc.c = Ac - 2Xb + TUc, (pcb,b = Ac - 2Xb + rzib, and (fcb = <fcb,b - <fcb,c = - ^c- For each 
of these angles, we measure the amplitude as half the difference between the maximum and 
minimum values obtained during each ~82 1-year simulation. As expected, we find that for 
both two- and three-planet fits to the real data set and for every bootstrapped data set, 
all three angles librate about 0°. We use the subscript "cb" to indicate the angles relating 
the longitudes of planets "c" and "b." In § [5l we will use the subscripts "be" and "ce" to 
indicate the corresponding angles relating the longitudes of planets "b" and "e," and "c" 
and "e," respectively. 

As in R05, we inclined the coplanar system to the line-of-sight. We set the inclination 
relative to the plane of the sky of all the companions to values from 90° to 40° in decrements 
of 1°. We hold all the longitudes of the ascending node fixed at 0°. We fit for all the 
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remaining parameters. Finally, we examine xt ^ function of the system's inclination. 
Figure [T] shows the results. 

The minimum in xt ^ function of the system's inclination is in the range 57-61°. 
We find that the location of the minimum has a small dependence on the number of 
companions in the system. For two planets, we find i = 57 ± 5.8°. For three planets, we 
find i = 59 ± 3.2°. For four planets, we find i = 61 ± 2.4°. (The results for four-planet 
fits are presented in § [5]). The error estimates on the coplanar inclination are derived from 
100 bootstrap trials of the RV data set. For each bootstrap RV set, we use the best-fit 
parameters to the real RVs at each inclination value as the initial guess. For each trial, 
we find the minimum in xt ^ function of the system's inclination. The uncertainties 
in the fitted inclination of the system are the standard deviations of the locations of the 
minimum in xt- Our best three- and four-planet dynamical fits indicate a system with an 
invariable plane with an inclination relative to the plane of the sky of ~ 59°, and we adopt 
this inclination as a starting point for studying the fitness of models with mutually inclined 
orbits. 

Starting from the coplanar, i = 59°, two-planet fit, we examine the effect of a mutual 
inclination between companions "b" and "c." We do this by fitting for the inclinations 
of "b" and "c" {ib and ic) and the longitude of the ascending node of "c" (fic)- Since 
the system is rotationally invariant about the line-of-sight, we leave the longitude of 
the ascending node of "b" [Q^) fixed at 0°. We find xt ~ 9-76 (rms=5.704 ms^^) for 
the coplanar fit, and xt ~ 9-91 (rms=5.706 ms~^) for the mutually inclined fit. This 
"improvement" is not significant. As we show below, fitting for the mutual inclination of 
"b" and "c" in the case of four-planet configurations results in a similar "improvement" 
in xt- For the two-planet mutually inclined model, the fitted masses and inclinations are 
rric = 0.7211 Mjup, = 2.3689 Mw, = 57.94°, and ib = 55.29°, and = 0.82°. With 
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the approximation that f2c ~ 0°, the invariable plane of this two-planet system has an 
inclination relative to the plane of the sky of 56°. The mutual inclination for this system is 
2.7°. 

For the three-planet model with i = 59°, our best coplanar fit has xt — 3-93 
(rms=3.625 ms~^), and the mutually inclined fit has xl = 3.99 (rms=3.620 ms~^). For 
the mutually inclined fit, the invariable plane of the system has an inclination relative to 
the plane of the sky of 59.2°. For the mutually inclined fit, we force planet "d" to be in 
the invariable plane determined by "b" and "c." The mutual inclination between "b" and 
"c" in this configuration is 0.7°. Since a mutual inclination does not significantly improve 
the three- planet fit, in Table [2] we show the best-fit parameters for our best coplanar 
three-planet fit. 

At z = 59°, the angles (pcb,c, ^cb,bi and (^cb are all in libration for both the model fit 
(in Table [2]) and for all of the bootstrap trials. For the fit in Table [2l these amplitudes 
are Lpcb,c = 4.5 ± 0.7°, (pcb,b = 13.1 ± 3.2°, and ipcb = 11-5 ± 3.7°. Note that the libration 
amplitudes for the three-planet fits presented in this section are smaller than for the 
three-planet fits presented in R05. 

The critical angles show a complicated evolution for the coplanar three-planet fit with 
i = 59°. The first angle, (pcb,c, has peaks at (listed with decreasing power) 1.58, 45.67, and 
8.69 years. The second angle, ipcb,b, has peaks at 8.69, 1.58, and 45.74 years. Finally, the 
third angle, (pcb, has two peaks at 8.69 and 45.67 years. These periodicities are also present 
in the evolution of the eccentricities of the two outer planets. The eccentricities of all three 
planets show very regular secular evolution. 
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Fig. 1. — For GJ 876, xl ^ a function of the coplanar system's inclination to the plane of 
the sky. The different types of points indicate how many planets are assumed to be in the 
system. Triangles are for two planets, squares are for three planets, and pentagons are for 
four planets. The eccentricity of "d" is allowed to float for all of the three- and four-planet 
fits used to make this figure. Note that the minimum is better constrained with the addition 
of each planet. 
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Residuals of the Three-Planet Fit 



We performed an error-weighted Lomb-Scargle periodogram ( IGilliland &: Baliunas 



19871 ) analysis on the residuals of the coplanar three-planet fits for both i = 90° and i = 59°. 
Both periodograms show a significant peak at ~125 days. Figure |5] shows the periodogram 
for the case i = 59°. The periodogram for the case i = 90° looks very similar. For i = 90°, 



the peak raw power is 311.4. For i = 59°, the peak raw power is 310. 2. Even if we a l 



Correia et al. 



ow for 



(120091) 



a mutual inclination between "b" and "c" and an overall trend as in 
(they actually subtract out a trend, due to the geometrical effect of the star's proper motion, 
prior to fitting the RVs), we still see a strong signal near 125 days in the periodogram of 
the residuals. For both values of i, folding the three-planet residuals at the peak period 
shows a strong coherent sinusoidal signature. Figure 13] shows the folded residuals for the 
case i = 59°. Examination of all the inclined three-planet fits shows that the peak period in 
the residuals is at 125 days for all values of i > 40°. For the 1000 fits used in determining 
the uncertainties in Table 121 the residuals of 941 of them have a peak period near 125 days, 
with 12 harboring a marginally higher alias peak with a period near 1.0 days. The 125-day 
periodicity has been adequately sampled at all phases only relatively recently. 

We also examined subsets of the RV data where we took just the first observations, 
assumed i = 59°, and fit for the parameters of the three known planets. We took all values 
of from 98 to 162. We then examined the error-weighted Lomb-Scargle periodogram of 
the residuals for all values of A^. All the periodograms have a peak period near 125 days. In 
Figure 13] we show the folded residuals for the first 98 observations as squares. Comparison 
between the folded residuals for the first 98 points and that for the full RV set again 
suggests that the 125-day periodicity has been sampled at all phases only recently. The 
observed power at 125 days has grown monotonically during the last decade, but the growth 
rate has been non-uniform. The power grows more quickly whenever we observe around the 
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time when the fourth planet is near the phases in its orbit when its orbital velocity has a 
large component along the line-of-sight. The power grows less quickly whenever we observe 
the system when the orbital motion of the fourth planet is mostly perpendicular to the 
line-of-sight. We are constrained to observe near the time of the full Moon, which impedes 
the efficiency of obtaining uniform phase coverage for periodicities near 120 days. As a 
consequence, there were several stretches where we consistently observed near the phases 
when the fourth planet's orbital motion was primarily perpendicular to the line-of-sight. 
This circumstance is one of the reasons why so many observations over 12 years were 
required prior to reliably detecting the fourth periodicity in this system. 

As an additional check on the reality of the 125-day signal, we ran a Monte Carlo 
analysis. We take the coplanar three-planet model at i = 59° and integrate forward to 
generate a synthetic RV curve. We then generate a set of 1000 synthetic data sets by 
repeatedly sampling the synthetic RV curve at the times of the Keck observations and 
adding Gaussian noise perturbations of amplitude equal to the rms of the fit in Table [2l 
Then, for each RV set, we fit for the three planets and examine the periodogram of the 
residuals. In not one instance did the power near 125 days exceed the value observed for 
the real data. Additionally, the peak period is "near" 125 days in only one of the 1000 
cases. Thus, the false alarm probability (FAP) is < 0.001. 

An alternative to the fourth planet hypothesis is that the 125-day periodicity is 
due to rotation of the star itself. The rotation period of GJ 876 is at least ~ 40 days, 
based on its narrow spectral lines and its low chromospheric emission at Ca II H&K 



( jPelfosse et al. 



19981 ). R05 presented photometric evidence of a rotation period of ~97 
days. Our measurements of the Mt. Wilson S activity index also confirms this rotation 
period. Thus, rotational modulation of surface features cannot explain the 125-day period 
in the velocities. The periodogram of our measured Mt. Wilson S activity index after 
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subtracting out a 14-year spot cycle has a peak period near 88 days. 



5. Four-Planet Fits 

The 125-day periodicity in the residuals to our best coplanar three-planet fit motivates 
exploration of the possibility that a fourth companion exists in the system. For the past 
several years, as we accumulated new RV measurements, we monitored the quality of 
dynamical fits containing a fourth planet with P ~ 125 days. From late 2005 through 
2008, these fits preferred a fourth planet with Cg > 0.2, and produced system configurations 
which, when integrated, were dynamically unstable. Only recently, have the four-planet fits 
resulted in systems which are stable on time scales exceeding a few hundred million years 
(at least for coplanar configurations). A fourth planet, with a mass similar to Uranus, and 
on a stable, low-eccentricity orbit, has now definitively emerged from the data. 

A primary reason why the latest fits result in a stable system is that we have only 
recently gathered enough observations to sample all phases of the 125-day periodicity. 
Figure [3] shows that most of the (small) phase gaps have now been filled in. Additionally, 
near the points with larger residuals, we also have additional observations with smaller 
residuals. As a result, the LM routine no longer prefers a fourth planet with period ~125 
days with a relatively large eccentricity. Another reason why it has been difficult to find 
stable four-planet fits with the fourth planet at 125 days is that the region around the outer 



2:1 MMR with companion "b" contains the boundary be tween stable and unstab 



for the two-planet system consisting of only "c" and "b" (iRivera fc Haghighipour 



e orbits 



20071). 



To highlight the generally tenuous stability of planets with periods near 125 days, we 
performed simulations based on the fit in Table [2] in which we added 1200 (massless) test 
particles on circular orbits in the plane of the system with periods in the range 120-125 
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Fig. 2. — Error-weighted periodogram of the residuals of the three-planet coplanar fit, with 
i = 59°, to the RV data of GJ 876. 
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0.5 1 

Orbital Phase 

Fig. 3. — Residuals of the three-planet coplanar fit, with i — 59°, to the RV data of GJ 876. 
The residuals have been folded at the period of the tallest peak in the periodogram of the 
residuals. The first 98 observations are shown as squares. 
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days every 0.5 days plus at the fitted period for "e." At each p eriod we placed 100 test 



partic les equally spaced in initial longitude. In agreement with 



Rivera fc Haghighipour 



( 120071 ) ■ nearly all are unstable on time scales < 10 years. Only one particle survived more 
than 10^ years; however, it was also lost after 7.4 Myr. It started with a period of 122.5 
days (astrocentric) at a mean anomaly of 219.6°. For this long-term survivor the critical 
angle, v^Lapiace = Ac — S\b + 2Ae, librates around 0° with an amplitude of A(y9Lapiaco = 22° 
for the entire time that the particle is stable. The chaotic nature of the particle's orbit 
becomes apparent when the libration amplitude of the Laplace angle suddenly becomes 
unconstrained during the last few thousand years of its evolution. In Figure S] we show the 
stability of the test particles as a function of initial semi-major axis and longitude. For 
clarity, we show the astrocentric semi-major axes raised to the sixth power. Unstable initial 
locations are indicated with open dark grey symbols. The single "stable" initial location, 
where a test particle survived 7.4 Myr, is indicated with a fiUed-in black symbol. With 
an open star, we also mark the approximate fitted location for the fourth planet from the 
coplanar four-planet fit with i = 59° (for simplicity, we assume zero eccentricity for the 
fourth planet — the astrocentric eccentricity is actually 0.0448). We also find that if we 
assign a mass of a few Earth masses to these "test particles," they tend to survive for longer 
times. The results strongly suggest that a resonant mechanism is required for an object in 
this region to be stable. 

Additionally, we find that the orbital parameters of the innermost planet, "d," can 
infiuence the stability and chaotic nature of test particles placed in the region near 125 
days. Assuming a circular orbit for "d" results in an island of stability near the fitted 
location of planet "e." Instead of only a single particle which survives for 7.4 Myr, six 
particles survive for at least 10 Myr in this region if is assumed to be zero. As for the 
long-lived particle for the case with > 0, all the "stable" particles are in the Laplace 
resonance. The six stable particles under the assumption = are shown as filled-in light 
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grey sym bols in Figure HI We also used the Lyapunov estimator in the SWIFT integration 



package (ILevison fc Duncan 



19941 ) to analyze the chaotic nature of test particle orbits with 
periods near 125 days. We find that the most stable particles in this region have Lyaponov 
times in the range 10^ to 10'^ years. 

Following procedures similar to the ones used for the two- and three-planet fits, we 
worked up four four-planet fits. The best-fit coplanar inclination {i = 59°) with = 
has xl = 2.8738 (rms=3.0655 ms~^). Working from this fit, holding fixed, but fitting 
for a mutual inclination between "b" and "c" gives xl = 2.8322 (rms=3.0219 ms~^). If 
instead we fit for and assume a coplanar system, the best-fit inclination is at z = 61° and 
xl = 2.5991 (rms=2.9492 ms^^). Finally, fitting for and a mutual inclination between 
"b" and "c" results in xl = 2.6098 (rms=2.9303 ms^^). Since fitting for and the mutual 
inclination between "b" and "c" results in a significant improvement in xl^ "we take this 
last model to represent our current best fit to the present Keck RV data set. However, 
simulating this model as well as other mutually inclined models with similar values of xl 
shows that the resulting systems are unstable on time scales of < 1 Myr. Also, with 
allowed to fioat, the "improvement" in xl between the coplanar and mutually inclined 
models is not significant. The F-Test probability that these two fits are statistically the 
same is 0.88. Thus, Table |3] shows the best-fit parameters for the coplanar four-planet fit 
with i = 59°, our preferred inclination for the system. The resulting system is stable for at 
least 300 Myr, and other coplanar fits with similar values of xl result in systems which are 
stable for hundreds of millions of years up to at least 1 Gyr. 

For completeness, we note that for our best-fit mutually inclined system, the invariable 
plane of the system is inclined 59.5° to the plane of the sky, and the mutual inclination 
between "b" and "c" is 3.7°. Fitting for the inclinations of planets "d" and "e" does not 
improve xl significantly. This result for the inclination of "d" is in agreement with the 
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similar result by ICorreia et al. I (120091 ) . The F-Test probability for our best four-planet 
(mutually inclined) fit versus the three-planet coplanar fit with i = 59° and with a fitted 
Cd is 0.0066. This small probability indicates that the four-planet model is significantly 
different from the three-planet model. It should be noted that for the coplanar case with 
i = 59° the fitted longitude of planet "e" is near the location of the long-lived particle 
discussed above. 

The inclination of planet "e" does have a significant effect on the stability of the 
system based on our best mutually inclined fit. It also influences the libration amplitudes of 
the critical angles to be discussed below. The orbit of planet "e" must have an inclination 
and node such that the precession rates of both the nodes and longitudes of periastron 
of the outer three planets result in a relatively small libration amplitude for the critical 
angle of the Laplace resonance, v^Lapiace = K — 3Afe + 2Ae. In other words, systems with 
a small forced inclination for "e" should be more stable than systems with a large forced 
inclination - through experimentation we find this to be true. The coplanar four-planet 
fits are generally more stable than the mutually inclined versions since a small libration 
amplitude for v^Lapiace depends only on the precession rates of the longitudes of periastron of 
the outer three planets and not on the rates of precession of their nodes. For the mutually 
inclined four-planet fit, we performed a brief search of the parameter space spanned by ie 
and Qe as well as and Q^- Experimentation shows that all four of these can influence the 
stability (and evolution) of this chaotic system. Also, the libration amplitudes of all the 
librating critical angles are smaller than if we simply place "e" and "d" in the invariable 
plane determined by "b" and "c." Since we only explored a small region of the parameter 
space, and since we found several equally good fits to the RV data with various values for 
the forced inclination of "e," other, more stable systems with smaller forced inclinations are 
likely to exist. One emerging pattern from the search of the parameter space is that smaller 
forced inclinations, smaller libration amplitudes, and more stable systems occur more likely 
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for models in which the orbital plane of "e" is more closely aligned with that of "b." This 
pattern is also in rough agreement with the finding that coplanar systems are generally 
more stable than mutually inclined ones. 

As for the two- and three-planet fits, the critical angles involving "c" and "b" librate 
about 0° for all our four-planet fits as well as for all 1000 bootstrap fits used in determining 
the uncertainties in Table [31 However, their amplitudes are generally larger because of the 
perturbing influence of the fourth planet. 

Since the ~125-day period is near the 2:1 MMR with planet "b," the 4:1 MMR 
with planet "c," and the Laplace resonance with both "b" and "c," it is informative to 
examine all the relevant critical angles associated with these resonances as well as those 
associated with the resonances involving only "b" and "c." For the 2:1 MMR and the 
linear secular coupling between "b" and "e," the relevant angles 

'^be,e = K — 2Ae + cCg, and (pbe = '^be,e — '^be,b = — For the 4:1 MMR between "c" 
and "e," there are 4 relevant angles: ipceo = K — 4Ae + SzUc, fcei = K — 4Ae + 2zuc + TUe, 
V'ce2 = Ac — 4Ae + + 2cCe, and (pce3 = K~ 4Ae + Sw^.. The subscripts above correspond to 
the multiplier in front of We- There is also, the critical angle involving just the longitudes 
of periastron of "c" and "e": ifce = ~ ^c- Finally, the critical angle for the Laplace 
resonance is v^Lapiace = Ac — 3Afe + 2Ae. Table H] lists the libration amplitudes of all these 
critical angles for the fit in Table [31 A "C" in the amplitude column indicates that the 
angle circulates (at least once) during the ~82 1-year simulation. The time evolution of the 
planetary eccentricities and the Laplace angle for the fit in Table [31 are shown in Figures \5\ 
and [6l respectively. The secular evolution of the eccentricities of planets "b," "c," and "d" 
is apparently slightly less regular as a result of the addition of planet "e," of which the 
eccentricity evolution is more chaotic than that of the other three planets. Additionally, the 
evolution of all the librating angles in stable four-planet fits appears chaotic; however, the 
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resulting systems are stable for at least several hundreds of millions of years. 

The addition of the low-mass planet "e" to the system model has only a modest 
effect on the dynamical characterization of the orbits of "b" and "c." A<^cf),c is increased 
from 4.5° ± 0.7° to 5.74° ± 0.85°, and A(/?efe increases from 11.5° ± 3.7° to 22.5° ± 4.8°. 
(Additionally, their evolution appears less regular with the addition of "e.") In addition, 
the 4:1 argument, v^ceo, is observed to librate with Ayjceo = 83° ± 25°, and the 2:1 argument 
ipbe,b librates with Ayjbe.fe = 36° ± 13°. Most significantly, the system also participates in the 
three-body resonance, with AyjLapiace = 40° ± 13°. 

Libration of v^Lapiace is required for stability. Long-term simulations based on bootstrap 
fits associated with four-planet coplanar fits indicate that systems for which (/^Laplace ^ 140° 
are unstable on time scales < 2 Myr. Other systems with large libration amplitudes for 
<^Lapiace, 141.4°, 129.1°, and 119.4°, are stable for only 25.9 Myr, 33.6 Myr, and 196.8 Myr, 
respectively. This suggests that the stability of the system also hinges on the libration 
amplitude of the critical angle for the Laplace resonance. 

Figure [7] charts the planetary positions over 120 one-day intervals starting from JD 
2450608.093 in a frame rotating at planet "b"'s precession rate, (wb) = —0.116165° d^^. 
Four snapshots, corresponding to successive half-orbits for planet "b" are shown. The plots 
trace the system's 1:2:4 commensurability at a time when the osculating eccentricity of 
"e," at Ce ~ 0.064, is relatively high, and when a/fc — ~ 180°. During this phase, triple 
conjunctions occur when planet "e" is near apastron. 

For the dynamical configuration given in Table [3l numerical integration indicates that 
the periastron longitudes zuc and tUf, regress more quickly than does We- Over a timescale 
of several years, the three apses approach alignment, angular momentum is transferred to 
"e," and the eccentricity of "e" 's orbit declines nearly to zero. During this low-eccentricity 
phase, the arguments ipbe and (pce both circulate through 2tt. The complicated dynamical 
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environment for "e" results in dramatic changes in its eccentricity over readily observable 
time scales. Our fit to the data indicates that during the time span of the RV observations, 
the eccentricity of "e" has varied from Cg ~ to Cg ~ 0.075 (see Figure E]). 

6. Are the Three Outer Planets Truly in a Laplace Configuration? 

In the previous section we showed evidence for a fourth planet in the GJ 876 system, 
planet "e," which participates in a 3-body resonance with planets "c" and "b." Our 
preferred coplanar model for the GJ 876 planetary system which is consistent with the RV 
data shows that the critical angle for the Laplace resonance librates about 0° with amplitude 
40°. Given the low RV half-amplitude for planet "e," it is natural to question whether there 
exists sufficient data to truly confirm the existence of the three-body resonance between 
"b," "c," and "e." In this section, we address the questions concerning the reality of this 
dynamical configuration. The analysis in this section is another method to address the 
question of whether the 125-day period is actually due to a planet. 

We performed a Monte Carlo simulation in which we took a three-planet model and 
added a ~ 125-day signal plus Gaussian noise. The three-planet model is just the three 
inner planets from the fit in Table [3l We added a Keplerian model with P = 120.713 days, 
m = 14.629 Me, e = 0.0448, u = 251.36°, and MA=216.88° plus Gaussian noise with 
rms=5.9696 ms~^. These (astrocentric) parameters are the result of subtracting off the 
effect of the three inner planets from the fit in Table [3] and fitting the residuals with a 
one-planet model (with initial guess based on the parameters of the fourth planet in our 
preferred four-planet fit) for which only the mean anomaly is allowed to fioat. We generated 
1000 synthetic RV data sets and performed a four-planet Newtonian fit to each data set. 
We then integrated each fit for 300,000 days (~821 years), and examined the libration 
amplitude of the critical angle for the Laplace resonance. Note (again) that the synthetic 
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a6 (AU) 

Fig. 4. — Stability of massless test particles placed at equally spaced initial longitudes on 
circular orbits with periods of 120-125 days in the plane of the system based on the fit from 
Table El For clarity, we show the astrocentric semi-major axes raised to the sixth power. 
Unstable locations are indicated with dark grey open symbols. The single "stable" location 
is indicated with a black filled-in symbol. The fitted initial location of the fourth planet from 
the fit in Table |3] is marked as an open star. Stable locations under the assumption 6^ = 
are indicated with light grey filled-in symbols. 
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Fig. 5. — Eccentricities of the four planets from the fit in Table |3] over the first 100 years of 
an 821-year simulation. 
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Fig. 6. — Critical argument for the Laplace resonance for planets c, b, and e over the first 
100 years of an 821-year simulation. The simulation is based on the parameters from Table O 
Horizontal lines indicate the observed amplitude for this simulation when it is extended to 
10^' years for x =5, 6, 7, and 8, respectively. The amplitude grows stochastically with time. 
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T = 60 days T = 90 days 

Fig. 7. — Four configuration snapshots of tlie GJ 876 planetary system. Eacli panel 
shows the positions of the four planets at 120 successive one-day intervals starting from 
JD=2450608.093 (T=0 days), with the orbital positions at the listed times given by red 
dots. The diagrams are drawn in a frame that rotates to match the mean orbital precession 
of planet "b," which amounts to —10.45° over the 90 days shown. Planet "b"'s apsidal line 
coincides with the x-axis. The apses for planets "c" and "e" are drawn with smaller-dashed 
and larger-dashed lines, respectively. 
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RVs are based on a Newtonian three-planet configuration plus a Keplerian signal plus 
Gaussian noise. Thus, we may expect that not too many fits should result in systems with 
a libration amplitude comparable to the system based on our preferred four-planet model. 
We find that only 6 simulations have a libration amplitude < 40°, and 29 have a libration 
amplitude < 53°, the one-sigma uncertainty on the Laplace angle libration width. Thus, 
most of these synthetic RVs result in systems which are not comparable to the system 
based on our preferred four-planet model. These results suggest that we can constrain the 
libration amplitude of the critical angle for the Laplace resonance. The results also indicate 
that planet "e" is strongly interacting with "b" and "c." This is in good accord with the 
eccentricity evolution shown in Figure |5l 

Another Monte Carlo experiment we performed is to assume the system is deep in the 
Laplace resonance and generate synthetic RVs to see how often we find systems comparable 
to the system based on our preferred four-planet fit. Among the fits to the bootstrap 
RVs above, we found one which results in a system in which the libration amplitude of 
the critical angle for the Laplace resonance is 16.4°. We take this as our model and add 
Gaussian noise with rms=3.9604 ms~^, the rms for the fit in Table [3l We generate 1000 
synthetic RV data sets and perform four-planet Newtonian fits to them. The 1000 fits 
are then integrated forward for ~ 821 years. We find that 577 systems have a libration 
amplitude of the critical angle for the Laplace resonance < 40°, 799 systems have this 
libration amplitude < 53°, and none have the critical angle in circulation. The largest 
amplitude is ~ 130°. Again, these results indicate that we are constraining the libration 
amplitude of the critical angle of the Laplace resonance. 

Additionally, we examined the long-term stability of several bootstrapped and synthetic 
systems for which the critical angle for the Laplace resonance is not strictly librating. We 
find that such systems are unstable on timescales < 200 Myr. Thus, our preferred coplanar 
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fit to the current data set shows a system which is participating in a Laplace resonance, 
and this configuration appears to be required for long-term stability. 

As noted in Figure |6l the libration amplitude for the critical angle for the Laplace 
resonance grows stochastically for long-term simulations performed up to at least 1 Gyr. If 
this diffusion occurs in the actual system and if we assume that the system were damped 
down initially, then the currently observed libration amplitude for the critical angle for the 
Laplace resonance may give a rough estimate of the age of the system. 



7. Are There More Bodies in the System? 



Rivera fc Haghighipour I ( 120071 ) used massless test particles to show that, except for 
the innermost region, the entire habitable zone is unstable for small planets because of 
the perturbing influence of just planets "c" and "b." Based on their work, it is expected 
that interior to the orbit of "e" there may be a small region of stability just exterior to 
the orbit of "d." Exterior to the orbit of "e," we expect the region out to the location of 
the 2:1 MMR with "e" to be unstable. Beyond this location, we expect small planets on 
low-eccentricity orbits to be stable except possibly around the location(s) of the 3:1 MMR 
(and perhaps the 4:1 MMR) with "e." 

We have performed a similar analysis based on the parameters in Table [3l In fact, due 
to the strong interactions among planets "e," "b," and "c," we flnd that particles placed 



exterior to the orbit of "d" are unstable out to b eyond the location o 



Rivera fc Lissauer 



the outer 5:2 MMR 



20001) in which the 



with "e." This result is similar to that found by 
interactions between the outer planets in the T Andromedae planetary system clear out 
a region much larger than what would be cleared out if only the outermost planet were 
present. 
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The analyses above were done with test particles placed on circular orbits starting 
on the a;-axis. They suggest that the region among the planets is unstable. This is also 
suggested by the inherently chaotic nature of the four-planet system. A more thorough 
exploration of the initial positions of test particles can be used to find amo ng the planets 



small islands of stability which are protected by some resonant mechanism. 



Correia et al. 



( I2OO9I ) postulated that a small planet could exist with a period near 15 days. Such a planet 
is not detectable with the current RV data set(s). The addition of the fourth planet, "e," 
may have a significant effect on the stability of such a planet. However, we explored the 
possibility that such a planet could be participating in a four-body resonance with the three 
outer planets. Short-term simulations of particles placed in this region do indeed show that 
a small fraction of them are participating in the putative four-body resonance. A few of 
these test particles have the relevant critical angle, = A/ — 4Ac + 5Ab — 2Ae, executing 
small amplitude libration, ~ 30°, about cf) ~ ±70°. This configuration would keep the 
potential fifth planet away from the line joining the outer three planets and the star when 
triple conjunctions occur among "c," "b," and "e." Again, the present RV data set(s) do 
not support such a planet, but it would be interesting to discover a resonant configuration 
which has no Solar System analogue. 



8. Summary and Discussion 

Using 162 RV observations taken at the Keck I telescope with HIRES, we have shown 
strong evidence for a fourth periodicity that is present in the radial velocity of GJ 876. This 
likely corresponds to a planet "e" with period ~124 days and a minimum mass of ~13.4 
Mq. This planet joins the previously known planets "d," "c," and "b" with periods of 
1.9378, 30.1, and 61.1 days, respectively. Assuming that the four-planet system is coplanar, 
we have shown that the inclination of the system relative to the plane of the sky is ~ 59°. 
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Our best self-consistent fit shows that the period and mass of the new companion are 
124.2 days and 15.4 M^. The value of above 1.0 for our fit in Table [3] suggests that 
the stellar jitter may be at the 1-2 ms^^ level and holds out the possibility of additional 
low-amplitude objects in the system which may be revealed with additional observations. 

We have also shown that this new planet likely participates in a complicated set of 
resonances involving planets "c" and/or "b." Simulations indicate that it is in a 2:1 MMR 
with planet "b," a 4:1 MMR with planet "c," and a Laplace resonance with both "c" 
and "b." Comparison with prior four-planet fits using fewer RV observations indicates 
that as more RVs have been accumulated, the "fitted" libration amplitudes of the critical 
angles associated with these resonances have generally decreased. Additionally, long-term 
simulations based on our fits indicate that the complicated dynamical structure is necessary 
for the system's long-term survival. It is important to note that our current best-fit coplanar 
model may change as we and other groups gather more RV data on this fascinating system. 
Thus, the system may be deeper in the resonances mentioned above. Additionally, future 
RVs may eventually constrain the mutual inclination between planets "b" and "c." 

Phenomenal signal-to-noise and a decade-plus observational baseline combine to make 
the GJ 876 system a touchstone for studies of planetary formation and evolution. The 
presence of well-characterized planet-planet resonances can potentially allow a detailed 
connection to be made between the properties of GJ 876's protoplanetary disk and the 
planetary system that eventually emerged. Clearly, the system was prolific, both in terms 
of the mass and the total number of planets produced. 

A zeroth-order prediction of the core-accretion paradigm for giant planet formation 
is that the frequency of readily detectable giant planets should increase with both 
increasing stellar metallicity and with increasing stellar mass (Laughlin et al. 2004, Ida 
& Lin 2004, 2005). During the past decade, both of these trends have been established 



- 33 - 



obs ervationally (see , e.g 



and 



Johnson et al. 



Fischer &: Valenti 



fl2008h and 



Bowler et ah 



( 20051) for a discussion of the metaUicity trend 



(120 lOl ) for discussions of the mass trend.) 



Until recently, however, there appeared to be little evidence for the strong expected 
planet-metallicity correlation among the ha ndful of M-dw a rf sta rs that are known to harbor 



giant planets. In particular, the results of 



Bonfils et al. 



(120051) suggested that GJ 876 



has subsolar metaUicity. This result naturally induces speculation that GJ 876's giant 
planets might be the outcome of gravitational instability (e.g. Boss 2000) rather than core 
accretion. 



Johnson fc Apps 



2009) 



developed by iBonfils et al. 



have recen tly provided an update to the metaUicity calibration 



(120051). The 



Johnson fc Apps 



(120091 ) calibration indicates that 



the planet-bearing M-dwarfs do appear to be systematically metal-rich, suggesting that 
there is no breakdown of the planet-metallicity correlation as one progresses into the red 
dwarf regime. In particular, the new calibration indicates [Fe/H] = 0.39; GJ 876 is more 
than twice as metal-rich as the Sun. 

A supersolar metaUicity for GJ 876 dovetails with a history of vigorously efficient 
planet formation, but there remain a number of fascinating questions with regards to 
the formation and migration mechanisms that permitted the final configuration to be 
assembled. We are preparing an in-depth analysis of these issues, even as we continue to 
gather more observations of this most eminently productive star. 
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Table 1. Measured Velocities for GJ 876 (Keck) 



BJD 


RV 


Unc. 


(-2450000) 


(ms-i) 


(ms-l) 


602.09311 


329.19 


1.79 


603.10836 


345.30 


1.81 


604.11807 


335.99 


1.92 


605.11010 


336.00 


1.93 


606.11129 


313.94 


1.88 



Note. — Table 1 is presented in 
its entirety in the electronic edition 
of the Astrophysical Journal. A por- 
tion is shown here for guidance re- 
garding its form and content. 
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Table 2. Three-Planet coplanar fit for GJ 876 with i = 59° 



Parameter 



Planet d 



Planet c 



Planet b 



P (days) 

a'^ (AU) 
K (ms-l) 
e 

a;(°) 
MA (°) 
offset (ms" 
Fit Epoch 

xl 

RMS (ms- 

'Pcb,b 
<Pcb 



(JD) 



1.937781 ± 0.000024 
6.78 ± 0.51 M® 
0.02080665 ± 0.00000018 
6.60 ± 0.47 
0.257 ±0.070 
229 ± 28 
358 ± 29 



30.0880 ± 0.0091 
0.7175 ± 0.0042 Mjup 
0.129590 ±0.000026 
88.72 ± 0.52 
0.25493 ± 0.00080 
48.67 ± 0.82 
295.0 ± 1.0 
50.7 ±0.4 
2450602.093 
3.9280 
3.6254 
4.5 ± 0.7° 
13.1 ± 3.2° 
11.5 ±3.7° 



61.1139 ± 0.0084 
2.2743 ± 0.0059 Mjup 
0.208301 ± 0.000019 
213.86 ±0.55 
0.0292 ± 0.0015 
50.7 ±4.5 
325.4 ±4.6 



''Quoted uncertainties in planetary masses and semi-major axes do not incorporate the uncer- 
tainty in the mass of the star 



Table 3. Four-Planet coplanar fit for GJ 876 with i = 59° 



Parameter Planet d Planet c Planet b Planet e 



P (days) 

a'' (AU) 
K (ms-l) 
e 

MA (°) 
offset (m s^-"-) 
Fit Epoch (JD) 

xl 

RMS (ms-i) 



1.937780 ±0.000020 
6.83 ±0.40 Me 
0.02080665 ± 0.00000015 
6.56 ± 0.37 
0.207 ± 0.055 
234 ± 20 
355 ± 19 



30.0881 ± 0.0082 
0.7142 ± 0.0039 Mjup 
0.129590 ± 0.000024 
88.34 ± 0.47 
0.25591 ± 0.00093 
48.76 ± 0.70 
294.59 ± 0.94 

51.06 ± 0.30 
2450602.093 
2.6177 
2.9604 



61.1166 ±0.0086 
2.2756 ± 0.0045 Mjup 
0.208317 ± 0.000020 
214.00 ± 0.42 
0.0324 ± 0.0013 
50.3 ±3.2 
325.7 ± 3.2 



124.26 ±0.70 
14.6 ± 1.7 Me 
0.3343 ± 0.0013 
3.42 ± 0.39 
0.055 ± 0.012 
239 ± 22 
335 ± 24 



''Quoted uncertainties in planetary masses and semi-major axes do not incorporate the uncertainty in the mass 
of the star 



Table 4. Libration Amplitudes of Relevant Critical Angles 



Angle Amplitude (°) 



<Pcb,c 


5.74 ± 0.85 


'Pcb,b 


21.9 ±4.2 


<Pcb 


22.5 ±4.8 


'Pbe,b 


36± 13 


Vbe,e 


C 


'Pbe 


C 


fceO 


83 ±25 


Vcel 


C 


'Pce2 


C 


VceS 


C 


fee 


C 


V'Laplace 


40± 13 
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